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1 INTRODUCTION 



The Zel'dovich 
IBuchert" 



(N 



approximation (ZA) (IZel'dovich 1970 
Goetz 19871 IBuchert 1992] IBouchet et al. 19951 
IBuchert fc Ehlers 19931 IBernardeau et al. 20021 

|Rampf fc Buchert 2012] |Rampf 2012| ) provides a very 
simple analytical model of the gravitational evolution 
of CDM inhomogeneities which reproduces the appear- 
ance of the cosmic web, correlating well with the large 

and void regions emerging in 



cosmic 

scale filamentary features 



non- linear N-body simulations ( Sathyaprakash et al 
|Springel 2005T 



IBuchert et al. 19971 

|Neyrinck 2012| |Rampf fc Wong 2012[ ) with the same initial 
conditions. The ZA can be derived from the full system of 
(Newtonian) gravitational equations and forms a subclass 
of solutions in Lagrangian perturbation theory (LPT) 
(Buchert 1992). In fact, the ZA and its second order 
improvement (2LPT) are used to provide the initial dis- 
placements and velocities of particles in N-body simulations 
(Scoccimarro 1998; Croccc, Pucblas & Scoccimarro 2006). 

The ZA arises in Newtonian theory and one might won- 
der about its status within general relativity (Elli s 20021 
IBuchert fc Oste rmann 2 012|l . This question is particularly 
relevant if the ZA is used for example to set the initial dy- 
namics of very large simulations which approach or exceed 
the size of the horizon. In this letter we show how the ZA 
and the 2LPT displacement field are derived in a general rel- 
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ABSTRACT 

We show how the Zel'dovich approximation and the second order displacement field of 
Lagrangian perturbation theory can be obtained from a general relativistic gradient 
expansion in ACDM cosmology. The displacement field arises as a result of a second 
order non-local coordinate transformation which brings the synchronous/comoving 
metric into a Newtonian form. We find that, with a small modification, the Zel'dovich 
approximation holds even on scales comparable to the horizon. The corresponding 
density perturbation is not related to the Newtonian potential via the usual Poisson 
equation but via a modified Hclmholtz equation. This is a consequence of causality 
not present in the Newtonian theory. The second order displacement field receives 
relativistic corrections that are subdominant on short scales but are comparable to 
the second order Newtonian result on scales approaching the horizon. The corrections 
are easy to include when setting up initial conditions in large N-body simulations. 

Key words: cosmology: theory - large scale structure of Universe - dark matter 



ativistic framework for ACDM cosmology. They correspond 
to a gradient expansion solution of the Einstein equations 
(Rigopoulos & Valkcnburg 2012]) , expressed in a coordinate 
system in which the metric takes a Newtonian form. In the 
process we calculate the relativistic corrections to the dis- 
placement field as well as the time shift between the proper 
time of the irrotational CDM particles and the "Newtonian" 
time corresponding to a weakly perturbed metric. 

As is the case in LPT, the gradient expansion allows 
in principle for density contrasts that are larger than unity, 
5p/p > 1. The expansion eventually breaks down at points 
where caustics occur and the density becomes infinite. Close 
to such singularities higher order terms in the gradient se- 
ries become important and the expansion loses its predictive 
power. However, one would expect that, unless black holes 
form, such singularities are in some sense removable since 
they appear where the worldlines of CDM particles cross. 
We find that when the gradient expansion breaks down, the 
corresponding Newtonian frame spacetime can still be con- 
sidered a weak perturbation of FLRW. 



1995 

IMatsubara 2008 



2 THE GRADIENT EXPANSION METRIC 

The gradient expansion is a technique for approximating so- 
lutions to the Einstein equations which is not based on ex- 
panding in small perturbations, as in conventional perturba- 
tion theory, but on writing the time-evolved metric in terms 
of a series in powers of the initial 3-curvature. The idea dates 
back to Lifshitz & Khalatnikov (1963) and Tomita (1975), 
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and was developed more recently in Stewart et al. (1994), 
Comer et al. (1994) and Tanaka & Sasaki (2006) (see also 
Comer (1996) and Bruni & Sopuerta (2003) for covariant for- 
mulations) . We will use here the gradient expansion solution 
for an irrotational flow of CDM particles in the presence of 
A (Rigopoulos & Valkcnburg 2012) to derive the Zel'dovich 
approximation^ 

Let us begin by writing the metric in synchronous co- 
moving coordinates, possible to construct in this case, 



As 2 



-At 2 + jij{t, q) dg l dg J 



(1) 



Summation over repeated spatial indices is implied. 
Here t is the proper time of the CDM particles and 
q are comoving coordinates, constant for each CDM 
fluid element. The metric can then be approximated by 
URigopoulos fc Valkenburg 2012) 



7»j(o,q) ^ akij + A(a) [ 



+ a 



Ax 



X(x)J{x) 
x 5 H(x) 



-a 2 I Ax- L{X) 



x 5 H(x) 



~rR kij 

4 



4R*. 



\2R kl R M k 



kl '» ij 



10R kl R k ik 



ISRRij + 32RikR 



a 

-a 2 j Ax 



X(x)J(x) - L(x) 
x B H(x) 



kij—4Rij' -k-\-R;i 



(2) 



In this expression kij is an initial "seed" conformal metric 
describing the geometry early in the matter era. We assume 
this initial conditions to hold as a — > 0, effectively setting the 
lower limit of the integrals in ([2]) to zero. Terms containing 
decaying modes have been neglected. Hats indicate that the 
curvature tensors are to be evaluated from the initial time- 
independent conformal metric kij, e.g. R = k tJ Rij(kti), and 
a semicolon "; k" denotes a covariant derivative w.r.t. this 
metric. We have used the background FLRW scale factor 
a(t) as the time variable and 



H(a) = H ^n M a-3 + Q A . (3) 
In terms of the proper time t of the CDM particles we have 

a(i) = expj^ At'H(t')\ , (4) 

with 



/3A 



H(t) = J — coth 
The functions appearing in the integrands in ([2]) satisfy 



(5) 



dJ J _ 1 
da a 2aH 



dA 
da 



, A 



J 

aH 



AL 
da 



L 
a 



J 



„ ■ (6) 
a atl da a an 

It is easy to write down the solutions to these as integral 
expressions but it is simpler numerically to solve the above 

1 See Barrow &: Goetz (1989) for a Newtonian treatment. 



equations directly. At early times, when the contribution 
from A is negligible, we have 




L 



,9/2 



175 H 3 D. 3/2 



(7) 



and a ~ t 2 ^ 3 (Hq-\/£Im)'' ' . These expressions are exact for 
an EdS universe with SI a/ — > 1 and Ho — > 2/(3to). 

Let us now focus on the following conformal seed metric 



kij — Sij 



l + y*(q) 



(8) 



where $(q) is the initial Newtonian potential, taken to be a 
Gaussian random field with amplitude given by the appro- 
priate transfer function. The metric kij is simply the linear 
initial condition derived from inflation and expressed in the 
synchronous gauge. Of course, non-Gaussian initial condi- 
tions could also be incorporated in kij but we keep (JSJ for 
simplicity. Solution (J2Jl becomes 



+ jA(o) 



1+— *(q) 



10 A 



+ Ti(a)$,ti<f> M ■ -T 2 (a) 



T 3 (a) 



where 



Ti = a' 



T 2 = a' 



and 



da: 200 
x 5 H(x)~ 

Ax 400 
x 5 H(x)~9~ 



X(x)J(x) - -L(x) 



X(x)J(x) - -L(x) 



(9) 

(10) 
(11) 

(12) 



A ", r denotes a differentiation w.r.t. Lagrangian coordinate 
qi. Note that a similar solution in the context of second 
order perturbation theory for Einstein-De Sitter cosmology 
was given in Matarrese, Pantano & Saez (1994). However, 
the expressions are not identical to ours; we have retained 
terms with two spatial gradients (the first two lines of ©) 
which are up to two powers in the potential $. These terms 
are crucial for the coordinate transformation below. 
The energy density is given by 

1 3/2 



p(a,q) = 



3i?o^ii 



+ f $( q ) 



8^G ^det[ 7lJ (a,q)] 



(13) 



which matches to the linear perturbation theory density in 
synchronous gauge at sufficiently early times. It should be 
stressed that in deriving ((2| no assumption has been made 
about the magnitude of the density perturbation and values 
of Sp/p > 1 are in principle allowed. Of course, the density is 
accurate only up to the gradient order kept in the expression 
for the metric. However, the metric in the form ([2jl or ((9) 
predicts that eventually regions of zero volume will form 
where the density becomes infinite^ As we will see below, 

2 See Rigopoulos & Valkenburg (2012) for the spherical case and 
the corresponding Zel'dovich approximation formula for the den- 
sity. 



© 0000 RAS, MNRAS 000, L1-L5 



Zel'dovich approximation and General Relativity L3 



this in general corresponds to the formation of caustics in 
LPT. 



3 NEWTONIAN COORDINATES AND THE 
DISPLACEMENT FIELD 

The spatial coordinate system used above is comoving with 
the CDM fluid, i.e., each fluid element, or particle, is charac- 
terised by a fixed q throughout the evolution. All informa- 
tion about inter-particle distances and clustering is encoded 
in the metric. This however is not the most convenient way 
to visualise the situation, to relate to Newtonian intuition 
or, for example, to compare with the output of an N-body 
simulation. Let us therefore define a coordinate transforma- 
tion from the comoving coordinates (t, q) to another coor- 
dinate system (r, x) where we require the metric to take the 
Newtonian form 



ff00 (r,x) =-[l + 2A(r,x)] , 






(14) 


goi(r,x) = 0, 






(15) 


ffij(T,x) =6ij [l-2B(r,x)]a 2 (r), 






(16) 


where A <C 1 and B <C 1, and where 








x\t,cfi = q +J 4 (t,q L ) , 






(17) 


r(t,ci)=t + C(t,q). 






(18) 


The metrics are related through 








dr dr „ ,, dx l dx m . 

™ = a* a/ 1 + 2A > + 


(i 


- 25) a 2 , 


(19) 


dr dr „ „. Ox 1 dx m r 


(i 


- 25) a 2 , 


(20) 


drdT „ n ^ dx l dx m . 

otm ( 1 + 2 ^ + at at Sl - 


(i 


- 25) a 2 . 


(21) 



In the above equations the various functions are evaluated 
at the same spacetime point, which we choose at this stage to 
label with the (t, q) coordinates that parametrize the world- 
lines of the CDM particles. For simplicity we will ignore pos- 
sible vector and tensor modes that are generated at next to 
leading order - this can be straightforwardly rectified. We 
can now obtain the displacement field and the time shift at 
different orders in the potentials 



F = F 1 (t,c i ) + F 2 (t,ci) + ... , 
£ = £i(t,q) +£ 2 (*,q) + ■■■ ■ 

3.1 Zel'dovich approximation 

Solving (I19|l - (|21[) at linear order in <E> we obtain 



(22) 
(23) 



m,q) = HJ||JU( q)j £ 1 (t, g ) = |j(t)$(q) 1 (24) 
and the gravitational potentials read 



Ai(t,x) = 5i(r,x) = ^[2H(t)J(t) - l]$(x) . 



(25) 



We see that the transformation (|17[) has a direct interpreta- 
tion: When expressed in terms of r it is simply the trajectory 
in the Newtonian frame (r, x) of a particle with initial co- 
ordinate q: 



. , 10 A(r) d , 

x(r, q) ~ q + — -±-L—$(q) , 
3 a 2 (T) 9q 



(26) 



where the replacement t — > r only induces a change at sec- 
ond order. We have checked that the prefactor 10A/(3a 2 ), 
although satisfying apparently different equations is numer- 
ically identical to the ACDM growth factor D+ (r) 



10 A(t) 
3 a 2 (r) 



D+(r) 



W(a) 



5n 



3/2 21 I 2' 6' 6 ' 



3 5 



o 
11 



dx 
W{x) 



A 3 

— a 



(27) 



with T-L the conformal Hubble parameter. The representa- 
tion of 5+ in terms of the hypergeometric function 2F1 was 
found in Enqvist & Rigopoulos (2011) and Belloso, Garcia- 
Bellido & Sapone (2011). We have thus obtained directly 
the Zel'dovich approximation for ACDM from a general rela- 
tivistic solution. Note that the formal steps described above 
resemble a gauge transformation from the synchronous to 
the Newtonian gauge. However, we stress again that we have 
not assumed here that Sp/p is smaller than unity. So, eq. (f26l) 
applies also in principle when <5p/p > 1. 

Focusing on scales that are comparable to the Hubble 
length, we can expand (| 13p to linear order in the potential 
and express the result in terms of r using (|24l) : 



-={r, x) 



A(t) 10 
p ' a 2 (r) 3 



V 2 $(x) + 10H(r) J(r) $(x) 



(28) 



We see that in the Newtonian frame and on scales compa- 
rable to the horizon the Newtonian potential and the den- 
sity perturbation are related through a (modified) Helmholtz 
equation instead of the standard Poisson equation. Writing 
the equation in terms of an evolving Newtonian potential 
<f>N (r, x) we have 

V 2 <Mt,x)-3^^<Mt,x) = ^, (29) 
A p 



with the non-local solution 



<Mt, x) 



dV 



exp 



2 HJi 



|x-y| 



5p 



47r|x - 



To be concrete let us set Qa = 0, Q n 
exactly. We then have 

3a 2 HJ 3H2 



A 



(r,y)- (30) 
■- 1 so that eq. J7J holds 



(31) 



We see that density fluctuations at distances sufficiently far 
away do not contribute to the potential. Furthermore, the 
region over which density fluctuations do contribute to the 
potential grows with time. This is of course expected since 
equation ([28} and the underlying solution ((26} are derived 
from general relativity. Such causal behaviour is absent in 
the Newtonian theory which misses the second term on the 
LHS of ([29}. 

Let us finally see how formula (|28p can be understood 
in terms of the particle trajectories (|26[l . Suppose for the 
moment that the Zel'dovich displacement (|26|) is imposed on 
a Euclidean grid. This would result in a density fluctuation 



s 4 

P / Euclidean 



A(r) 10 
a 2 (r) 3 



V 2 $(x). 



(32) 



The true spatial geometry of the Newtonian frame is not 
Euclidean and the density contrast acquires an extra term 
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(10HJ — 5) $ due to the change of spatial volume associ- 
ated with (|25[1 . Comparing with (|28f) we see that to obtain 
the correct density a condition on the initial Lagrangian po- 
sitions of the particles must be imposed. Indeed, assuming 
particles initially displaced by Xi n i = q + c(q) with 

V, ■ c = -5$ . (33) 

and evolved with (|26p will reproduce the correct density. 
This is in agreement with the result of Chisari & Zaldarriaga 

(2011) . A detailed discussion about the initial conditions 
can be found in Buchert (2011), Buchert, Nayet & Wiegand 

(2012) and Rampf & Rigopoulos (2012). 



3.2 Trajectory at second order and short-scale 
behaviour 

The displacement field, the time shift and the gravitational 
potentials can be calculated to second order as well. Explicit 
expressions can be found in the appendix. The trajectory in 
the Newtonian frame (r, x) of a particle with initial coordi- 
nate q reads at second order 



x(r,q) 



10 X(t) d . . . 
3 a (r) dq 

50 [\(t) + J 2 (t)] d 



+ l_T 2 (r) 8 



■ a 2 (r) dq V 



9 a 2 (r) 
50 [2A(r) + J 2 (r 



<9q V 



3 _t_ 
2 V 2 , 



F 



9 



a 2 (r) 



(iq 



(34) 



where 1/ V, denotes the inverse Laplacian. The last term 
in the first line of (|34[1 is precisely the result from 
Newtonian 2LPT (Rampf fcBuchert 2012) |Rampf 2012) 
|Rampf fc R igopoulos 2012 I. On small scales the second line 
in (|34[) is completely negligible and we obtain the com- 
plete Newtonian result, showing that Newtonian dynam- 
ics on short scales produce the correct evolution. This is of 
course not surprising, since general relativity is constructed 
to reproduce Newtonian physics in the appropriate limit. 
However, on scales approaching the horizon the last two 
terms in (|34[1 become comparable to the second order New- 
tonian terms, the ratio between the two scaling roughly as 

Relativists . . «| Thig ghQwg ^ firgt order m _ 



Newtonian k' 2 

suit, at second order general relativistic effects do have an 
impact on the trajectories of particles on such scales. In 
particular, any deviation from the Zel'dovich approximation 
computed with Newtonian dynamics on scales approaching 
the horizon will introduce errors. However, since dynamics 
on such scales are accurately described by expression (|34|l , it 
is easy to include these corrections in an N-body simulation. 
We give more quantitative details in Rampf & Rigopoulos 
(2012). 

Finally, let us make a few comments for the regime 
where the gradient expansion breaks down. On short scales 
the second order potentials read (see appendix) 



A 2 (r,x) = B 2 (t,x 
25 1 



9 a 2 1 



+ 



(2XHJ — A — J 2 — HL) + C($ 



(35) 



where we have dropped terms that are not enhanced by spa- 
tial gradients; denotes differentiation w.r.t. Eulerian co- 
ordinates xi, and F is the analogue to F in Eulerian space. 



It is interesting to examine what happens to the metric po- 
tentials when the gradient expansion solution breaks down. 
To simplify the expressions let us set Qa = 0, fi m = 1 so 
that eq. holds exactly. Expression ([9]) predicts that the 
components of the synchronous metric will go to zero for 
fluctuations with highest wave-number k approximately at 
a time defined by 



4 k 2 
3H? 



(36) 



At this point the spatial volume element of the comoving 
synchronous hypersurfaces goes to zero and the density be- 
comes infinite. In the Newtonian frame this signifies the 
crossing of the CDM particle worldlines and the formation 
of caustics (shell crossings). At these points the gradient 
expansion solution breaks down. 

Ultimately, shell crossings are the result of the assump- 
tion of a single velocity for each fluid element. In reality such 
singularities will be smoothed out by the non-zero velocity 
dispersion of the CDM particles but the zero pressure gradi- 
ent expansion solution used here will be inaccurate in these 
regions. However, some qualitative estimates can be made. 
At the time when these singularities form we approximately 
have 



A 2 



Bo 



84 



(37) 



We thus see that the second order correction to the metric is 
enhanced, formally becoming first order. This would signify 
that, even if zero pressure is still assumed, the complete se- 
ries should be summed to obtain the correct spacetime met- 
ric. However, unless terms of successive orders become even 
more dominant, eq. (|37|) shows no evidence that spacetime in 
these high density regions will be significantly different from 
FRW. This statement of course would be incorrect close to 
the formation of black holes which cannot be seen in this 
formalism. But this should not be the case for most such 
regions. 



4 SUMMARY AND DISCUSSION 

We have shown that the application to our universe of the 
gradient expansion method for approximating solutions to 
the Einstein equations is the relativistic equivalent of solving 
Lagrangian Perturbation Theory. At first order the relativis- 
tic displacement field coincides with the Zel'dovich approxi- 
mation up to an extra initial displacement c(q) which has to 
be imposed on the initial positions of particles to reproduce 
the correct density. We have therefore found that even for 
scales close to (or larger than) the horizon, the Zel'dovich ap- 
proximation is essentially correct as a description of particle 
motion. However, the relation between the resulting density 
contrast and the Newtonian potential is not the standard 
Poisson equation but a modified Helmholtz equation, reflect- 
ing the causality of the relativistic theory. Contrary to what 
happens at first order, the second order displacement field 
receives relativistic corrections that are as important as the 
corresponding Newtonian result on large scales. 

One can draw two main conclusions from the above 
findings. The first is that the fully relativistic solution re- 
produces the Newtonian dynamics on short scales. This is 
of course not surprising. However, we believe this is the first 
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time it is explicitly shown starting from a fully relativis- 
tic solution, making no assumptions on the magnitude of 
the density perturbation. Our results show no evidence that 
Newtonian cosmology is not a good description on short 
scales. Correspondingly we expect any backreaction to be 
small even when large density contrasts form. This finding 
should be compared with the backreaction estimated in the 
synchronous gauge (Kolb et al. 2006; Enqvist ct al. 2012). 

The second conclusion is that on large enough scales rel- 
ativistic effects start contaminating the second order New- 
tonian result with the relative importance of the relativistic 
terms scaling as Ho/k 2 . Such corrections will be relevant for 
simulations that encompass the horizon. Since the Zel'dovich 
term will dominate on such scales, the corrections will be 
rather small. However, any deviation from the Zel'dovich ap- 
proximation computed purely through Newtonian dynamics 
will miss the relativistic corrections in (|34[) . Formula (|34[) 
then provides a direct way to include relativistic effects on 
the trajectory of particles in large N-body simulations. We 
will return to this issue with a more quantitative treatment 
in Rampf & Rigopoulos (2012). 
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APPENDIX 

In this appendix we summarise our findings of the second 
order transformation to the Newtonian frame. It reads 



100 A d ^ 2 
8a 2 d<f T Vl r 9 a 2 «"« 



50 , t2 ,s 1 d 1 
+ — {J+ A) 



3 1 



9 v " a? dq* V 2 V 2Vf 



(38) 
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+f(;-^)(i^*A-(^ 

-f (J + 2HJ 2 )* 2 , 
while the second order potentials are 
a , \ 25 x ^ 25/3 18 rT T 2 t2 > , ., 

+ f±(2XHJ~X-J 2 -HL)^F 
_^Q_3HJ +7 // 2 J 2 -J 2 A 



2 1 * 1 

X I o 



(39) 



F , (40) 



25 A 



25 2 



Ba(r ' x) = \T^ 1 ^ 1 + ~ < ^HJ-SH 2 J 2 + -J 2 A ) $ 



9 V5 



2 



+ (2XHJ — A — J 2 — HL) -^-F 

9 a 2 y ' V 2 

+ ^.HJ (1-2H J) ( % -^.^ - _L_F ] ill) 



3 V 2 



Note again that the dependences and derivatives in eqs. (|38[1 
and (|39|1 are w.r.t. Lagrangian coordinates, while for 
eqs. (|40|l and (|41|) w.r.t. Eulerian coordinates. 



© 0000 RAS, MNRAS 000, L1-L5 



